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I.   INTRODUCTION 

Let  Au  =  s  be  a  linear  system  of  algebraic  equations  resulting  from 
the  discretization  of  an  elliptic  boundary  value  problem.   Direct  methods  of 
solution  that  factor  A  into  a  lower  triangular  matrix,  L,  and  an  upper  tri- 
angular matrix,  U,  such  that  A  =  LU,  are  computationally  inefficient  since 
the  factorization  yields  unique  L  and  U  which  are  not  sparse  matrices.   For 
this  reason,  an  iterative  procedure  is  more  generally  desirable. 

Iterative  procedures  have  recently  been  proposed  by  several  people, 
notably  H.  L.  Stone  [3]  and  [k-~\,    and  Dupont,  Kendall  and  Rachford  [1],  based 
on  the  idea  of  forming  a  matrix  (A+B  )  which  can  be  factored  into  the  product 
of  sparse  triangular  matrices,  L  and  U  ,  i.e., 


(1.1)  A+B  =  L  U 

a    a  a 


The  solution  of  Au  =  s  is  then  computed  as  the  limit  of  the  sequence  defined 

by 

(1.2)  (A+Ba}Vl  =  (A+VVT(AVs)  • 

where  a  and  i   are  iteration  parameters.   The  procedures  of  Stone,  and  Dupont, 
Kendall  and  Rachford  differ  only  in  the  method  used  to  determine  (A+B  ) . 

Another  efficient  iterative  technique  for  solving  the  linear  system 
Au  =  s  is  the  alternating  direction  implicit  (ADl)  procedure.  ADI  was  ini- 
tially proposed  by  Peaceman  and  Rachford  [2].   The  idea  is  to  partition  A  as 
the  sum  of  two  easily  invertible  matrices  H  and  V,  i.e., 

(1.3)  A  =  H  +  V  . 


The  following  iteration  scheme  is  then  used  to  obtain  the  solution,  u: 

(H+wihl)un+l/2  =  S  "  (V-Wihl)un 

(V+u.  I)u    =  s  -  (H-w.  I)u  .  /r> 
iv   u+1  iv  '    n+1/2 

where  to.,  and  w.   are  iteration  parameters, 
lh      iv 

In  his  paper  [3],  Stone  described  a  factorization  procedure  that 
proved  to  have  several  advantages  over  those  now  in  use.  His  studies  indi- 
cated that,  in  addition  to  reducing  computational  effort,  the  convergence 
rate  did  not  depend  strongly  on  the  nature  of  the  coefficient  matrix,  A,  or  the 
shape  of  the  region  being  considered  and  suitable  iteration  parameters  could 
easily  be  calculated.   For  these  reasons,  his  method  is  the  subject  of  much 
current  research,  but,  so  far,  attempts  at  rigorous  analysis  have  been  unsuc- 
cessful.  One  of  the  difficulties  is  that  his  factorization  replaces  A  with 
an  asymmetric  matrix,  (A+B  ),  and,  moreover,  this  matrix,  being  dependent  on 
the  variable  parameter  a,  changes  frequently  during  iteration.  A  new  factor- 
ization recently  devised  by  Stone  [k]   defines  a  symmetric,  positive  definite 
matrix,  (A+B,.),  to  replace  the  asymmetric  one  in  the  procedure.   Since  sym- 
metry greatly  simplifies  matrix  algebra,  it  is  reasonable  to  expect  progress 
in  the  mathematical  analysis  of  this  symmetric  procedure.   If  so,  these 
findings  might  also  be  applicable  to  the  asymmetric  procedure.  However,  the 
symmetric  procedure  is  not  just  an  approach  to  the  asymmetric  procedure.   If 
it  is  competitive  to  ADI,  it  too  is  of  practical  interest.   Therefore,  the 
purpose  of  this  thesis  is  to  test  the  effectiveness  of  Stone's  symmetric  fac- 
torization by  comparing  the  results  obtained  with  those  from  both  Stone's 
asymmetric  procedure  and  the  ADI  procedure .   In  addition,  since  proper  usage 
of  parameters  is  important  to  convergence  rates,  pertinent  data  will  be 
compiled  to  help  determine  their  properties. 


II.   MATHEMATICAL  PRELIMINARIES 


Consider  the  Dirichlet  problem 

<«•«       -  4  (*!«%)  -  4  (a2(x%) =  s(x) ' xeD 

u(x)  =  f(x),  xeD 

where   =  (x  ,x  ),  D  is  the  interior  of  a  compact  region,  and  the  differential 
operator, 

<2-ia>  -  4  iai<x^)  -  4  va2(x)4) ' 

is  elliptic  in  D. 

In  order  to  solve  this  problem  numerically  we  will  let  D  be  the 
interior  of  the  unit  square,  see  Figure  1,  i.e., 

(2.2)  D={(X;L,xx)  :0<x1,x2<1], 

and  cover  it  with  n  horizontal  and  n  vertical  grid  lines  with  equal  spacing 
h  =  — =-  ,  n  >  0  and  let  D  ,  the  grid  system  over  D  with  uniform  spacing,  be 
defined  by 

(2.3)  I>h  b  ((ih,jh)  :  1  <  i,  j  <  n)  ;  i,  j  integral. 

The  grid  point  (in, jh)  will  be  denoted  by  (i,j)  and  the  boundary  dD  of  EL  by 


(h,  0) 


(nh,  0)    (1,0) 


Figure  1.      Grid  Point   System 


'(i,0)   ,   0  <  i  <  n  +  1 
(1,  j)   ,   0  <  j  <  n  +  1 


(2.4)  ^\  =    < 


(i,l)   ,   0  <  i  <  n  +  1 

(0,5)      ,   0  <  j  <  n  +  1 
s. 

We  will  replace  the  differential  operator  (2.1a)  with  the  difference  operator 

(2.7)  and  replace  u  with  a  vector  whose  components  correspond  to  the  grid 

points  of  the  system.   This  yields  a  set  of  simultaneous  linear  equations 

having  a  one  to  one  correspondence  with  the  grid  points. 

The  grid  points  will  be  ordered  in  two  ways  and  the  ordering  may 

alternate  from  one  iteration  to  the  next.   The  first,  called  normal  ordering, 

arranges  the  grid  points  in  a  left  to  right,  down  to  up  pattern.   Thus,  grid 

point  (i, j)  precedes  grid  point  (k, j)  if  i  <  k,  and  precedes  grid  point  (l,m) 

if  j  <  m.  The  second,  called  reverse  ordering,  arranges  the  grid  points  in  a 

left  to  right,  up  to  down  pattern.   In  this  way,  grid  point  (i,j)  precedes 

grid  point  (k, j)  if  i  <  k,  and  precedes  grid  point  (l,m)  if  j  >  m.   Thus, 

normal  ordering  is  down  to  up  and  reverse  ordering  is  up  to  down.  Note  that 

components  of  the  iteration  vector,  u.  .,  must  conform  to  the  ordering  of  the 

i,  J 

grid  points  at  all  times. 

If  u  is  the  iteration  vector  with  components  u.  .,  (i, j)  a  grid 

i,  J 

point,    let  A  defined  by 

(2.5)  (Au).     .   =  H.     .u.     .    n    +  F.     .u,    _     .   +  E.     .u.     . 

i,J  i,0   i,J-l  i,J   i-l,  J  i,J    i,J 


+  F.    n     ,U.  ..     .    +  H.     .    ..u.     .    . 
l+l,  J    l+l,  j  i,J+l   1,  J+l 


be  a  discretization  of  the  differential  operator  (2.1a).   The  important 


properties  of  A  are  that 


(i)  H   ,F    <  0  for  1  <  i,  j  <n  ; 


(ii)  F    =  0  for  i  =  1,  j  =  1,  ...,n  ; 

i,  0 

(iii)   H    =  0  for  j  =  1,  i  =  1,  . .  ,,n   ; 

i,  0 

(2.6)            (iv)   E.  .  =  -(H.  .  +  F.  .  +  F.  _  .  +  H.  .  ,)  >  0  for 

1,0      1,0  1,0    1+1,  G    1,0+1 

1  <  i,  j  <  n,   if  F    ^  0  and  H.  .  J  0   ; 

(v)   E.  .  >  -(H.  .  +  F.  .  +  F.  ,  .  +  H.  .  .,)   for 

1,0     i,j  i,o   i+i,o   1,0+1 


1  <  i.  j  <  n.   if  F.  .  =  0  or  H.  .  =  0, 
-  '       -    '  1,0  1,0 


Note:   If  i  =  n,  F.  ,.  .  =  0  and  if  j  =  n,  H.  .  ,  =  0 
'  1+1,0  '  1,0+1 


The  above  properties  are  a  result  of  many  possible  discretizations  of  (2.1a), 
for  example, 

(2.7)       -D_Xi(a1(x+e1  j&D^uW)  -  D^fa^x^  |)D^u(x)) 

where  e.  is  the  unit  vector  in  the  direction  of  the  positive  i   coordinate 
i 

axis  and 


(2.8)  hD   u(x)  =  ±(u(x±he.)  -  u(x)) 

ix .  l 

l 


This  discretization  yields 


H,    ,   =   -a.(ih,  (j-  hh)    ,    if     j   =  1,    H,     .   =  0      ; 


i,J 


i>J 


Hi,j+1=   -a2(ih,(J+|)h)    ,    if     j   =  n,    H.^.+1  =  0      ; 


F. 


.      1 


=   -a   ((i-  ±)h,jh)    ,    if     i  =  1,    F         .  0 


.      1 


i+l,J   =   "a!((i+2)h^h)    «    if     i  =  n>    Fi+l,d  =  °      ; 


E.     .    =    -(H.     .    +   F.     .    +   F.    _     .    +  H.     .    J       . 
i>J  1*J  i,J  1+1,  J  i,  J+l 


III.   DESCRIPTION  OF  TEST  PROBLEMS 

Four  distinct  problems,  identical  to  those  described  by  Stone  in 
[3],  wiH  be  studied.   Each  problem  will  be  solved  using  the  square  grid 
system  described  in  Figure  1.   The  solution,  u(x),  of  (1.2)  and  (l.^)  was 
chosen  to  be 

(3.1)  u.  .  =  sin(ifrh)  •  sin(  jTfh)   for  i,  J  =  1,  ...,n 

1,  J 

where  n  is  the  grid  size  and  h  =  -> — =-y.   The  data  vectors,  s,  of  (1.2)  and 
(lA)  were  computed  from 

(3-2)  (AuE).  .  =  s.  .  . 

1,0    1,0 

The  initial  vector,  u.  .,  was  chosen  equal  to  zero. 

1,0 

The  first  test  problem  is  the  model  problem,  defined  by  letting  the 
coefficients,  a  (x)  and  a  (x),  be  equal  to  -.5  over  the  entire  region  described 
in  (2.3). 

For  the  second  problem,  the  generalized  model  problem,  the  ratio  of 
a  (x)  to  a  (x)  is  equal  to  100,  e.g.  a  (x)  =  -.5,  a  (x)  =  -.005. 

The  third  problem  is  defined  by  a  heterogeneous  region  containing 
homogeneous  subregions;  see  Figure  2.   In  region  A,  the  ratio  of  a^(x)  to 
a  (x)  is  unity  with  a  (x)  =  -.5.   In  region  B,  this  ratio  is  ten  with  a-[(x)  = 
-.5;  in  C  it  is  one-tenth  with  a  (x)  =  -.05,  and  in  D  it  is  one  with  a-^x)  = 

-.05. 

The  last  problem  is  identical  to  the  third  with  one  difference.   In 
region  A,  the  values  of  the  a. (x)  were  randomly  chosen  from  the  set  {z:0  <  z  <  1 
subject  to  the  condition  that 


(0,  1)   30 


(0,0) 


25       _ 


20       _ 


*      15 

-a 


10       _ 


5       — 


Figure  2.       Heterogeneous   Test  Problem  Geometry 
(Region  Defined  on  Unit  Square) 


(3.3)  If     z  <  -1     set  a-(x)  =  10~ 


otherwise    set     ai(x)    =   z  ,      i  =  1,2 


10 


11 

IV.   THE  STONE  PROCEDURE 

Direct  methods  to  solve  the  linear  system,  Au  =  s,  involving  the 
L-U  decomposition  of  the  coefficient  matrix,  A,  are  computationally  inefficient 
when  the  grid  size  is  large;  because  L  and  U  will  possess  n  non-zero  diagonals. 
To  make  L-U  decomposition  efficient,  A  is  replaced  by  a  matrix  of  the  form, 
A  +  B  ,  which  can  then  be  factored  as  the  product  of  sparse  lower  and  upper 
triangular  matrices,  L  and  U  ,  respectively.  An  approximate  solution  of 
Au  =  s  is  then  computed  by  stopping  the  iteration, 

(4.1)  (A+B  )u  .  =  (A+B  )u  -  t(Au  -s)  , 

v  cr  n+1      an      n    ' 

provided,  of  course,  that  convergence  is  sufficiently  rapid  to  be  practical. 

In  order  to  increase  convergence  rates  it  is  plausible  that  the 
vector  norm, 

C*.2)  l|Bau||, 

should  be  minimized.   Stone  suggested  constructing  L  and  U  ,  see  Figures 
3  and  4,  to  be  sparse  matrices  with  three  non-zero  diagonals,  each  of  which 
coincide  with  the  non-zero  diagonals  of  A,  i.e., 


(*.3) 


(L  u) .  .  =  b.  .u.  .  _  +  c.  .u.  _  .  +  d.  .u.  . 
a  1,  j    i,  j  l,  j-1    i,j  l-l,  J    i,j  i,o 


(U  u) .  .  =  u.  .  +  e.  .u.  ,  ,  +  f .  .u.  .  , 
a  i,  o    i,0    i>J  i+l,  J    i,0  i,0+l  • 


Figure  5  and  equations  (4.3)  show  that  the  components  of  (L  U  u)  are  given  by 


12 


Figure  3.   Lower  Triangular  Matrix,  L     Figure  h.      Upper  Triangular  Matrix,  U 


a 


a 


Figure   5.      Product  Matrix  L^Uq,  =  A+^, 
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(L  U  )u  .  .  =  b.  .u.  ,  ,  +  b .  .e.  .  ,u.  ,  .  _  +  c.  .u.  ,  . 

-  ddf  Ji,j    1,0  i,o-l    i,j  1,0-1  1+1,0-1    i,  0  i-l,  0 

n.   \,\  +  (d.  .+b.  .f.  .  _+c.  ,e.  ,  .)u.  .  +  d.  .e.  .u.  .  . 

(4.4)  i,J   1,0  1,0-1   1,0  1-1,0   1,0     i,J  i,J  1+1,0 

+  c.  .f.  ,  ,vl.   .  ...  +  d.  .f.  ,u.  .  , 

1,0  1-1,0  i-i,  0+1    1,0  1,0  1,0+1 

where  the  subscripts  i, j  in  (4.3)  and  (4.4)  are  the  grid  point  indices,  not 
row-column  indices.   (Throughout  this  section,  lower  case  lettering  will  be 
used  to  represent  components  of  the  product  matrices,  L  and  U  .   Capitals 
are  reserved  for  labeling  the  components  of  the  coefficient  matrix,  A,  and 
the  auxiliary  matrix,  B  .)   Figure  9  displays  the  components  of  the  iteration 
vector,  u,  that  appear  in  (4.4)  and  the  following  equations  are  Taylor  series 
expansions  for  these  components  in  terms  of  the  value  of  u  at  the  grid  point 
(i,0): 

(i)  u  ,  ,.,  =  u.    -  hu „(i,o)  +  hu  (i,o)  +  0(h  ) 

p 
(ii)  u       =  u    -  hu  (i,o)  +  0(h  ) 

i-i,  0      x}  J     -*■ 

2 
(iii)  u.  .,.   su.  .  +  hu  (i,o)  +  0(h  ) 

1,0+1      i,0  y 


(4.5) 


(iv)   u.  .     =  u.  . 
1,0        1,0 


2 
(v)  u       =  u    +  hu  (i,j)  -  hu  (i,o)  +  0(h  ) 

i+i,o-i   i,  o    x       y 


(vi)  u       =  u    +  hu  (i,o)  +  0(h  ) 


2 

(vii)   u.  .  _    =  u.  .  -  hu  (i,o)  +  0(h  ) 

1,0-1      i,0  y 


(viii)   u.  .     =  u.  . 
1,0        i,0 


The  matrix,  B  ,  is  then  constructed  such  that:   the  elements  of  any  row  mul- 
tiply only  those  components  of  u  which  appear  in  (k.k)    and  (U.5);  and  the 
vector  norm,  (h .2) ,    is  minimized.   Figure  7  depicts  an  asymmetric  matrix,  B  , 


Ik 


u 


whose  components,  (B  u),  are  given  by 


(B  u).  .  =  -qC.  .u.  .  .  +  C.  .u.  _  .  _  -  aG.  .u.  ,  . 

a  i,j     1,0  1,0-1   i,  j  i+i,  j -l    1,0  1-1,0 


(+.6) 


+  a(C.  .+G.  .)u.  .  -  aC.  .u.  ,  .  +  G.  .u.  _  .  _ 
'1,0  i,  0  i,J     1,0  1+1,0    1,0  1-1,0+1 


-  aG.  .u.  .  _  . 
i,0  1,0+1 

Note  that  when  a  =   1.0,  the  quantity,  ±G.  .,  is  a  multiplier  for  the  first 

1,  J 

four  equations  in  (^.5)  and  the  quantity,  ±C.  .,  is  a  multiplier  for  the  last 

i,  0 

four  equations  in  (4.5)-   Consequently,  the  value  of  (^.2)  is  reduced  to  second 

2 
order,  0(h  ),  terms,  when  the  u.  .  are  given  as  values  of  a  smooth  function. 

i,  0 

The  algorithm  for  determining  the  elements  of  L  and  U  for  this 
to  to  a     a 

asymmetric  factorization  is  derived  by  equating  like  terms  of  the  matrices 
(L  U  )  and  (A+B  ),  Figures  5  and  8,  respectively.  The  resulting  equations 
are 

b.  ,  =  H.  ,  ,/(l+ae.  .   ) 
i,j    i,j-l     1,0-1 


c.  .  =  F,  -  ,/(l+af.  .  .) 
i,J    i-l, J      1-1,0 


(U.7) 


d.  ,  =  E.  .  +  a(C.  .  +  G.  .)  -  b.   f.  .    -  c.   e 
i,j    1,0      1,0    1,0     1,0  1,0-1    i,0  i-i,, 


e.  ,  =  (F.  .  -  aC.  .)/d.  . 

i,j      1,0      1,0    1,0 


f .  .  =  (H.  ,  -  aG.  . )/d.  . 
1,0      i,0      1,0    i,0 


15 


;'-QC»  C...-OG'    Cg£,  -QC...G'  -QG1 


Figure  6.   Coefficient  Matrix,  A 


Figure  7«   Auxiliary  Matrix,  B 


(Asymmetric) 


a 


.•       H.     .      -QC»      C...F 


-OG'      E.     .-fQC'+oG'      F.     .-QC...G1      H.     .-QG' 
1_1>J  i*J  i,J  i,J 


Figure  8.        Asymmetric  Altered  Matrix,   A+B 


Q 
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3>l,j 


i»J-l 


ifl,d-l 


Figure  9.      Points  of  Grid  System  Corresponding  to  the  j-th  Row  of 
the  Auxiliary  Matrix,    B 


a 
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where  i,j  =  l,...,n  and 


(U.8) 


C.  ,  =  b.  ,e. 


G.    =  c.   f    , 
i,j    i,j  l-l, j 

A  second  method  of  minimizing  (k.2)    is  obtained  by  constructing  the 

matrix,  B  ,  as  shovn  in  Figure  11.   In  this  case,  if  q.=  1,  (B  u)  is  first 

order  if  the  components,  u.  .,  are  given  as  values  of  a  smooth  and  continuous 

i  »0 
function.   This  factorization  has  the  improved  property  that  B   is  symmetric. 

The  following  equations  for  determining  the  elements  of  L  and  U  are  obtained 
&  ^  &  a     a 

by  equating  like  terms  of  the  matrices,  (L  U  )  in  Figure  5  and  (A+B  )  in  Figure 
12,: 


b.  .  =  d.  .  _f.  .  . 
i,j    i,J-l  1,0-1 


c.  .  =  F.  ,  .  -  aG  .  .  , 
1,0    i-l,0      1,0-1 


(^•9)  d.  .  =  E.  .  +  2aG.  .  .  -  b.  .f.  ,  _  -  c.  ,e,  .  . 

i,0    1,0       1,0-1    isO  i,0-1    i,0  1-1,0 


i.     .  =  (F.  .  -  aG.  .  .  _)/d. 
i,0     1,0      1+1,0-1   i, 


f.  .  =  (H.  .  -  aG.  .)/d.  . 
1,0    '1,0      i,0    1,0 


where  i,  j  =  1,  .  .  .,n  and 


G.  .  .  =  b.  .  .e 
1,3-1    1-1,0  1-1,0-1 


(+.10)        G.  ,  .  _  =  b.  ,e. 

l+l,  j -1    i,0  1,0-1 


G.  .  =  c.  .f.  .  . 
1,0    1,0  1-1,0 
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Figure  10.   Coefficient  Matrix,  A         Figure  11.   Auxiliary  Matrix,  B 


(Symmetric) 


Q 


H.  .   -0G.   G^...F.  ,  .-QGn   E.  .+2QGn   F   .-QG_...Gn   H.  .-QGn 
i,j-l   1   2    i-l,  J    1   i,J     1   i,j    2    3    i,j    3 


Figure  12.   Symmetric  Altered  Matrix,  A+B 
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If  u  is  not  smooth,  cancellation  is  not  effective  and  Stone  suggests 
using  a  variable  a  .   If  0  <  a  <  1,  then  whenever  the  matrix  modification 
introduces  the  value  u.  ,  .  .,  into  the  individual  grid  point  equation,  it  is 
to  be  partially  cancelled  by  adding  the  value  a(u.  .  -  u.     -  u    .).   Stone 
suggests  using  a  set  of  m  parameters,  a.,  i  =  1,  ...,m,    0  <  a.  <  1  computed  by 
letting: 


(4.11)    a  = 


m 


a 

max 


1  -  min 


2Axn 


a2  (x)  A  x£ 


1  + 


a1(x)Ax2  , 


2Ax 


ax(x)  Ax2 


1  + 


a2(x)  Axx 


and  o:  =  a  .   =0;  the  intermediate  parameters  a.,  i  =  2,  .  ..,m  -  1  are  computed 
from 


(4.12) 


i-1 

a.   =  1  -  (1  -  a   )m_1  ,  i  =  2, . .  .,m  -  1 

i  max     '  '        ' 


which  insures  that  0  <  a.    <  a.   .  <  a 


1,  . .  . ,  m  -  1 


i    "i+1  —  "max, 
For  the  test  problems  in  this  thesis  involving  Stone's  procedure 

with  parameters,  a  set  of  nine  parameters  was  computed  and  each  was  used  on 
two  successive  iterations.   For  all  odd  iterations  beginning  with  the  first, 
the  solution  was  obtained  using  the  normal  ordering  of  the  grid  points  as 
explained  in  Chapter  II.   For  all  even  iterations  beginning  with  the  second, 
the  solution  was  obtained  using  the  reverse  ordering  of  the  grid  points,  also 
explained  in  Chapter  II.  Altering  the  ordering  for  each  iteration  requires 
that  the  coefficients  of  L  and  U  ,  given  by  equations  (^.7)  and  (4.9),  be 
recomputed  each  time.   Since  each  parameter  was  used  on  two  successive  itera- 
tions, the  parameter  cycle  length  is  18.   Before  the  parameters  were  applied 
to  the  iteration  scheme,  they  were  numbered  in  order  of  increasing  magnitude, 
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i.e.,    cl    <  ql    <  . . .   <  a  ,    and  then  rearranged  in  the   following  manner: 

°ty  %>  Oy  Vq>  <y  <*2>  ^(\>a1- 

Stone  [3]  remarked  that  if  the  values  of  the  coefficients  E,  F,  and  H  in  the 
matrix  ,A,  see  Figure  10,  represent  sharp,  point  to  point  discontinuities,  and 
if  the  parameter  a  is  varied  during  the  iteration  scheme,  then  any  values  of 
x  other  than  unity  do  not  generally  improve  the  method.   Since  no  numerical 
or  analytical  justification  for  this  statement  appears  in  [3],  tests  were 
performed  using  other  values  of  x.   These  tests  show  an  improvement  and  the 
results  are  discussed  in  more  detail  in  Chapter  VI. 

Numerical  studies  to  determine  the  properties  of  the  iteration 
parameters,  a  and  x,  were  performed  using  the  symmetric  factorization .pro- 
cedure. For  each  test,  a   and  x  remained  fixed  for  a  maximum  of  32  iterations. 
For  each  change  of  a  or  x,  the  initial  vector,  u.  .,  was  chosen  as  zero.  A 
systematic  search  for  the  optimum  value  of  x,  i.e.,  the  value  that  provided 
optimal  convergence  for  32  iterations,  was  made  for  a  =   1.0,  in  the  range 

C  2.   For  the  model  problem  on  a  30  x  30  grid,  divergence  resulted  after 
a  few  iterations  whenever  a  was  less  than  unity.   Since  no  value  of  x  could 
effect  convergence,  a  was  fixed  at  1.0  for  the  remainder  of  the  tests. 
Figures  13  and  Ik   show  the  results  of  tests  performed  on  the  model  problem 
and  the  generalized  model  problem  for  various  grid  sizes  ranging  from  5x5 
to  30  x  30.   Similar  results  do  not  appear  for  problems  3  and  k   because  of 
prohibitive  work  and  time  requirements.   Only  the  optimum  values  for  x  are 
shown.  Although  only  grid  sizes  which  are  multiples  of  5  were  used,  a  con- 
tinuous curve  was  drawn  to  indicate  a  trend.   Figures  13  and  Ik   indicate  that 
(i)   the  optimum  value  of  x  decreases  as  the  grid  size  increases; 
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Figure  13 •  Model  Problem  -  Tau  vs.  Grid  Size 


1.0, 
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Figure  ik.    Generalized  Model  Problem  -  Tau  vs.  Grid  Size 
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(ii)  the  optimum  value  of  t  is  less  affected  by  the  grid  size  when 
the  ratio  of  a  (x)  to  a  (x)  is  large.   This  ratio  was  equal 
to  100  in  Figure  lk. 
For  all  problems  the  calculated  optimum  value  of  x  was  less  than  unity, 
regardless  of  grid  size. 

Other  tests  were  made  to  determine  the  effect  of  reordering  the 
grid  points  on  alternate  iterations.   Results  were  conclusive  in  showing  that 
more  effective  cancellation  and,  therefore,  faster  convergence  is  achieved. 
More  details  are  contained  in  Chapter  VI. 

The  actual  programming  of  the  procedure  described  in  (U.l)  was  done 
as  follows : 

Define 


and 


Solve 


R  =  t(s  -  Au  ) 

n  n 


6  ,  ,  =  u    -  u   . 
n+ 1    n+ 1    n 


L  V  =  R 
a  n 


for  V  and 


U  6  4.1  =  v 

a   n+1 

for  5  J_, .   Then  set 
n+1 

n+1    n    n+1  . 
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V.   THE  ADI  PROCEDURE 

A  well  known,  efficient  iterative  method  for  solving  large  systems 
of  algebraic  equations  of  the  form  Au  =  s  on  rectangular  grid  regions,  such  as 
described  in  (2.2),  is  the  alternating  direction  implicit  (ADI)  procedure.   The 
coefficient  matrix,  A,  is  partitioned  into  the  sum  of  two,  real  and  symmetric 
matrices,  H  and  V,  where  H  and  V  are  discretizations  of 


(5-D  -^V1'^ 


and 


(5.2)  -3x7(a2U)^) 


•> 


respectively.      The   solution,   u,   of  (H+V)u  =  s   is   computed  using  the   following 
2-step  iteration  scheme: 

(H+Vl)Vl/2  =   -(V-V)Un+   S 


(V+oj.  I)u  ..  =  -(H-oj.  I)u  .-  /0  +  s 
jv   n+1        jv   n+1/2 


where  u^  is  an  initial  estimate  of  the  exact  solution,  u,  u„  and  gj  .   are 
0  jh      jv 

iteration  parameters,  and  j  =  n(modulus  t),  where  t  is  the  parameter  cycle 
length.   In  (5»3),  2*t  iteration  parameters  are  used  in  a  cyclic  fashion  and 
changed  at  each  step  of  the  iteration  procedure. 

The  problem  we  are  concerned  with  is  how  to  choose  parameters  to 
minimize  the  computer  time  required  to  obtain  the  unknown  vector,  u,  to  some 
preassigned  accuracy.   The  solution  to  this  problem,  see  Wachspress  [6], 
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requires   that  the  matrices,   H  and  V,    commute  and  that  A  "be  positive  definite. 
If  these  two   conditions   are   satisfied,    optimum  parameters   that  minimize  the 
spectral  radius,   u(T , ),    of  the  t-step   iteration  matrix, 

*  -1  -1 

Tx   =      n    (V+w.    I)      (H-o).    l)(H+a)     T)       (V-oo.,1)       ; 
t         ,=1  jv  jv  jh  jh 


(5.M 


(ut_uE)    =  Tt(uQ-uE)       , 


can  be  easily  and  efficiently  computed  by  a  well-known  algorithm  [5],  which  cal- 
culates 2*t  optimum  parameters  by  minimizing  the  following  continuous  function: 


t   X-w.  Y-w., 

I       n  jv        jh 

y     =  max  n  u —  •   u— 


t  '     .    .X+w...  Y+w. 

J=l        Jh  jv 


(5.5)  a<X<b 


c    <_  y  £  d 

a  +■    c   =    X    .    (A)    >    0 
mm 


where  a,  b,  c,  and  d  are  eigenvalue  bounds  for  H  and  V,  respectively.   Since 
every  eigenvalue  of  T   is  of  the  form 


t   X-u.    Y-w.i 

(5.6)  n  — ^ & 

j=l    Jh     jv 


it   follows  that   after  t   iterations,   the   spectral  radius,    y(T    ),    is   such  that 


u,  -u 


(5.7)  y(T   )  =  Hi" 


<  u. 


V  ii  E.  |      -  Pt      * 

VU      1 2 
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No  scheme  has  yet  been  devised  for  choosing  optimum  parameters  when 
H  and  V  do  not  commute.  A  theorem  in  [2]  proves  that,  if  the  above  calculated 
parameters  are  used,  despite  lack  of  commutation  of  H  and  V,  any  desired 
error  reduction  can  be  obtained  if  the  cycle  length,  t,  is  sufficiently  long. 
It  is  our  experience,  see  also  Wachspress  [6],  that  repeated  application  of  a 
shorter  cycle  has  proved  better  when  H  and  V  do  not  commute.   Convergence  rates 
improved,  dramatically,  when  t  was  reduced  from  32  to  9,    while  the  number  of 
iterations  remained  32. 

Since  the  calculation  of  iteration  parameters  requires  the  values 
of  the  eigenvalue  bounds  for  H  and  V,  a  method  must  be  available  to  compute 
them.   The  values  of  the  upper  eigenvalue  bounds,  b  and  d  of  H  and  V,  respec- 
tively, are  less  critical  and  can  be  obtained  by  Gerschgoren' s  Theorem. 
Values  for  the  lower  eigenvalue  bounds,  a  and  c  of  H  and  V,  respectively,  are 
very  critical  for  good  parameter  calculation  and  must  be  determined  more 

carefully.   Wachspress  [7]  suggests  an  algorithm  based  on  Newton's  method 
applied  to  the  characteristic  polynomials,  |H-\I |  and  |V-\l|,  which  can  be 

evaluated  at  any  point,  X,  by  a  simple  recursion  [J] .  If  A  is  positive  defi- 
nite, it  follows  that  H  and  V  are  positive  definite  ,and  a  reasonable  starting 
value  for  Newton's  method  is  zero. 
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VI  RESULTS  AND  CONCLUSIONS 

Figures  15-l8  display  results  corresponding  to  the  four  problems 
described  in  Chapter  III  that  were  obtained  by  solving  each  problem  on  a 
30  x  30  grid  system.   In  each  of  the  figures,  the  ordinate  is  the  L  relative 
error  defined  by 


(6.1) 


where  u  is  the  exact  solution,  u^  is  the  initial  estimate,  and  u  is  the 

'     0  '     n 

estimate.   For  the  AD I  procedure,  two  half  steps  are  equivalent  to  one 
Stone  iteration.   The  abscissa  is  the  number  of  iterations  needed  to  obtain 
the  given  L  norm  ratio.   Only  straight  line  approximations  to  the  curves 
are  shown  in  order  to  facilitate  visual  comparisons.   The  actual  curves  con- 
tain oscillations  caused  by  parameter  cycles. 

Each  of  the  four  Figures,  15-l8,  contains  six  graphs  labeled  A-F. 
Each  letter  designates  one  of  the  following  methods  of  solution: 

A.  ADI  procedure  with  9  parameter  pairs  used  cyclically. 

B.  Stone's  asymmetric  procedure  with  a  parameters  used  cyclically 
as  described  in  Chapter  IV,  with  alternating  normal  and  re- 
verse ordering.  The  parameter,  t,  is  fixed  at  unity. 

C.  Same  as  B,  except  t  is  fixed  initially  at  the  value  that  pro- 
vides optimal  convergence  for  32  iterations.  In  all  figures, 
except  Figure  16,  t  is  1.6.   In  Figure  16,  x  equals  1.2. 

D.  Same  as  B,  but  with  Stone's  symmetric  procedure. 

E.  Stone's  symmetric  procedure  with  normal  ordering  on  every  iter- 
ation.  The  parameter,  a,  is  fixed  and  equal  to  unity.   The 


parameter,  x,  is  fixed  and  chosen  as  described  in  C  above.     27 
For  Figures  15,  16,  IT  and  18,  the  optimum  values  for  x  are 
.365,  .866,  .3^1  and  .326,  respectively. 
F.   Same  as  E,  except  with  alternate  normal  and  reverse  ordering. 
The  optimum  values  for  x  are  also  different.   They  are  .331, 
.78U,  .319  and  .32U,  respectively. 
By  studying  the  results  shown  in  Figures  15,  16,  17  and  18,  we 
are  able  to  make  the  following  observations  and  conclusions: 

(i)   For  all  test  problems,  procedure  B  converged  faster  than  pro- 
cedure D.   Therefore,  the  symmetric  factorization  is  not  as 
effective  as  the  asymmetric  factorization  of  Stone. 

(ii)   For  the  generalized  model  problem,  the  ACT  procedure  provided 
faster  convergence  than  any  of  Stone's  methods.   This  con- 
tradicts the  results  given  by  Stone  in  [3] •   The  reason  for 
this  discrepancy  is  that  Stone  did  not  use  optimum  iteration 
parameters.   Improvement  in  the  performance  of  the  ADI  pro- 
cedure is  not  as  dramatic  for  the  remaining  problems  as  for 
problem  2 . 
(iii)   In  every  case  procedure  C  out-performed  procedure  B,  disputing 
Stone's  statement  that  x  =  1.0  gives  better  results.   It  must 
be  noted,  however,  that  the  optimal  values  for  x  used  in 
procedure  C  were  found  by  a  tedious  systematic  search.   To 
fix  x  at  unity  is  as  reasonable  as  an  inefficient  search  for 
a  better  value . 

(iv)   The  results  of  procedure  F  are,  for  all  cases,  an  improvement 
over  those  of  procedure  E.   Since  the  only  difference  between 
these  two  procedures  is  the  ordering  of  the  grid  points,  the 
suggestion  by  Stone  [3]  to  use  alternate  normal  and  reverse 
ordering  to  increase  convergence  rates  is  numerically  justified. 
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(v)   Stone  [3],  in  an  attempt  to  show  that  one  of  the  sequence  of 
values  of  a  should  be  less  than  unity,  argued  for  large  grid 

sizes,  that  if  both  components,  F  and  H,  of  the  coefficient 

7rAx        IT ^  V 
matrix,  A,  are  large  relative  to  — ^ —   and  — — —     ,  respec- 
tively, then  repeated  application  of  a  =  1.0  in  the  iteration 
scheme,  (k.l),    leads  to  divergence.   Stone,  however,  fixes  x 
to  be  unity.   For  a  =  1.0,  values  of  x  other  than  unity  can 
be  selected  to  yield  convergence,  as  procedure  F  in  Figures 
15,  17 >    and  18  show, 
(vi)   The  results  of  procedure  F,  when  compared  with  those  of  pro- 
cedure D,  show  that  the  iteration  scheme  defined  by  (^.l) 
converges  faster  when  a  sequence  of  values  of  a  is  used.   The 
results  from  procedure  F,  for  problems  1,  3>  and  K,    are  optimal 
using  fixed  a  and  x  in  (^+.1)  because  experiments  indicate  that 
the  optimal  values  of  x  and  a  occur  when  a  =  1.0.  Values  of 
a  less  than  1.0  cause  divergence,  independently  of  the  value 
of  x  used  for  the  run.   Therefore,  finding  the  optimal  values 
of  (x,  a)  reduces  to  finding  the  optimal  value  of  (x,  1.0). 
Although  the  findings  in  observation  (i),  showing  the  symmetric  procedure,  D, 
to  be  less  effective  than  the  asymmetric  procedure,  B,  are  somewhat  discour- 
aging, they  do  not  conclusively  show  that  further  consideration  of  procedure 
D  is  useless.   It  is  more  important  that  the  results  of  procedure  D  be  similar 
to  those  of  B  as  the  following  observations  point  out: 

(vii)   For  the  complex  problems,  3  and  k-   in  Figures  17  and  l8,  proce- 
dure D  out-performed  the  ADI  procedure.   The  improvement  be- 
comes more  pronounced  as  the  complexity  of  the  problem  increases, 
see  Figure  l8.  Hence,  Stone's  symmetric  procedure  merits 

further  analytical  attention. 
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(viii)   For  the  simpler  problems,  1  and  2,  in  which  there  were  no  sharp, 
point  to  point  discontinuities  in  the  values  of  the  coefficients, 
F  and  H,  in  the  matrix,  A,  the  symmetric  procedure,  D,  was 
nearly  as  effective  as  the  asymmetric  routine,  B,  see  Figures 
15  and  l6.  Thus,  for  simple  problems,  any  findings  resulting 
from  analytical  studies  of  Stone's  symmetric  procedure  may  be 
applied  to  his  asymmetric  procedure. 
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